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Abstraot 


It  it  shown  in  this  paper  that  certain  orthogonal  polynomials  over  two 

disjoint  intervals  can  be  particularly  nsefnl  for  solving  large 

symmetric  indefinite  linear  systems  or  for  finding  a  few  interior 

eigenvalues  of  a  large  symmetric  matrix.  There  are  several  advantages 

of  the  proposed  approach  over  the  techniques  which  are  based  upon  the 

polynomials  having  the  least  uniform  norm  in  two  intervals.  While  a 

theoretical  comparison  will  show  that  the  norms  of  the  minimal 

polynomial  of  degree  n  in  the  least  squares  sense  differs  from  the 

minimax  polynomial  of  the  same  degree  by  a  factor  not  exceeding 
/l/2 

2(n+lK ,  the  least  squares  polynomials  are  by  far  easier  to  compute 

and  to  use  thanks  to  their  three  term  recurrence  relation.  A  number 
of  suggestions  will  be  made  for  the  problem  of  estimating  the  optimal 
parameters  and  several  numerical  experiments  will  be  reported. 
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l.  IfltrgdMtigB 

Th la  paper  i a  concerned  with  the  eolation  of  large  indefinite  linear  systems 

A  *  «  f.  The  noaber  of  effective  iterative  aethoda  for  treating  each  probleaa  ia 

2 

United.  An  obviooe  poeaibility  ia  to  aolve  the  noraal  equation*  A  x  «*  A  f  by  any 
available  aethod  for  aolving  positive  definite  aysteas  bat  this  is  not  reoooaended 
in  general  as  it  increases  the  aaoont  of  work  and  above  all  squares  the  original 
condition  noaber.  Perhaps  the  best  known  contribution  in  the  field  was  nade  by 
Paige  and  Saonders  who  adapted  the  Conjugate  Gradient  aethod  to  indefinite 
systeas  [8].  There  have  also  been  a  few  atteapts  towards  the  generalisation  of 
the  Chebyshev  iteration  algoritha  (also  oalled  Stiefel  Iteration)  to  the 
indefinite  ease  by  Lebedev  [6],Roloff  [10]  and  de  Boor  and  Bice  [4].  These  works 
focus  on  the  difficult  problea  of  finding  a  polynomial  of  degree  n  such  that 
Pn(0)~l  and  having  least  uniform  norm  in  two  intervals  containing  the  eigenvalues 
of  A  but  not  the  origin.  Like  in  the  positive  definite  case,  such  polynomials  are 
the  best  for  making  the  residual  norm  small  in  a  certain  sense  in  the  Richardson 
iteration: 

^  "I 

x  ^,-x  +x  r 
n+1  n  n  n 

where  t  is  the  n1^  root  of  the  polynomial  and  r  is  the  residual  of  x  .  There 
n  s&  n 

are,  however,  soae  serious  drawbacks  of  this  approach  one  of  the  aost  important 
being  that  the  aini-aax  polynomial  is  rather  difficult  to  obtain  nuaerically. 
Moreover  the  above  Richardson  iteration  is  unstable  in  general. 

At  this  point  one  aight  ask  whether  it  is  essential  to  use  the  ainiaax 
polynomial.  Perhaps  the  aost  attractive  property  of  the  Chebyshev  polynomials  is 
its  three  tera  recurrence  foraula  which  enables  one  to  generate,  in  the  poaitive 
definite  oase,  a  sequence  of  approxiaate  solutions  to  the  linear  system  which 
satisfy  a  three  tera  recurrence.  As  pointed  out  by  de  Boor  and  Rice  [4],  it  seeas 
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unfortunately  impossible  to  exhibit  t  three  ten  recurrence  for  the  minimax 
polynomial a  over  two  interval*. 

Nevertheless  a  natural  alternative  conaidered  in  thi*  paper  i*  to  construct  a 
sequence  of  polynomials  satisfying  p  (0)*1  and  which  minimise  an  L.  non 

A  2 

associated  with  a  certain  weight  funotion  over  the  two  intervals.  Such  polynomials 
are  easy  to  compute  and  can  be  obtained  via  a  sequence  of  polynomials  satisfying  a 

three  ten  recurrence.  Moreover  it  will  be  shown  that  their  unifon  non  in  the 

1/2 

two  intervals  is  comparable  within  a  factor  of  2(n+l)  with  the  non  of  the 
minimax  polynomial  if  n  is  the  degree  of  both  polynomials.  We  will  call 
Generalized  Chebyshev  Iteration  (GCI)  the  iterative  method  which  at  each  step 
yields  a  residual  vector  of  the  fon  r  -p  (A)r„  where  p  is  the  minimal 
polynomial.  There  are  several  other  applications  of  the  minimal  polynomial  in 
addition  to  the  Generalized  Chebyshev  Iteration.  Plrat.  a  block  version  of  the 
algorithm  can  be  fonulated.  It  is  essentially  a  generalisation  of  the  block 
Stiefel  algorithm  described  in  [11]  and  is  particularly  suitable  for  parallel 
computation.  Secondly  there  are  various  ways  of  using  the  minimal  polynomials  for 
computing  eigenvectors  associated  with  interior  eigenvalues.  Por  example  one  can 
easily  adapt  the  subspace  iteration  method  (see  [9])  to  interior  eigenvalues.  A 
few  of  the  possible  generalizations  and  extensions  of  the  basic  Generalized 
Chebyshev  Iteration  will  be  briefly  outlined  but  we  will  not  attempt  to  describe 
them  in  detail. 

Section  2  describes  polynomial  iteration  for  solving  linear  systems.  In 
section  3  we  will  deal  with  the  orthogonal  polynomials  over  two  intervals  and  will 
introduce  the  minimal  polynomial.  Section  4  covers  the  application  of  the  minimal 
polynomials  to  the  solution  of  indefinite  systems  while  section  5  outlines  their 
application  to  the  interior  eigenvalue  problems.  In  section  6  a  number  of 
praotioal  details,  eoneerning  in  partieular  the  estimation  of  the  optimal 


BP 
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parameter*,  are  given.  Finally  some  numerical  experiments  are  reported  in  section 
7  and  a  tentative  conolnsion  is  drawn  in  the  last  section. 


2.  Polynomial  iteration  is*  linear 


Consider  the  N  x  N  linear  syatem 


A  x  -  f 


(2.1) 


where  A  is  a  nonsingular  symmetric  matrix.  Let  Xq  be  any  initial  guess  of  the 

solution  x  and  rQ  the  corresponding  initial  residual  r0  ■  f  -  A  xQ.  Polynomial 

iteration  methods  are  methods  that  provide  approximate  solutions  x_  whose 

n 

residuals  r  satisfy: 

& 

rn  "  pn  (  A  >  r0  l2'2) 

where  p  is  a  polynomial  of  degree  n,  often  called  the  residual  polynomial.  We 
n 

only  consider  methods  in  which  x&  is  given  by 


n+1 


xn  +  °n°n 


(2.3) 


where  u  ,  the  direction  of  search*  is  a  linear  combination  of  the  previous 
n 

residuals.  In  this  case  it  can  easily  be  shown  by  induction  that  the  polynomial  p 

n 

satisfies  the  condition  p  (0)  **  1. 

n 


If  xQ  is  written  in  the  eigenbasis  (u^)^  N  of  A  as 


N 

*o  “  l  Vi 

i-i 


then  obviously  the  Euclidean  norm  of  r&  is 

1 1',"  ■  <  «i2  >,/2 


(2.4) 


where  the  X  ’  s  are  the  eigenvalues  of  A, 
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Equation  2.4  suggests  the  general  principle  that  if  we  want  the  residual  non 

Mr  ||  to  be  snail,  we  nnst  find  a  polynomial  satisfying  p  (0)  =  1  which  is  snail 
n  n 

in  sons  set  containing  the  eigenvalues  of  A.  One  way  of  choosing  a  residual 
polynomial  which  is  small  in  a  set  containing  the  spectrum  is  to  take  among  the 
polynomials  satisfying  p&  (0)  =  1  the  one  which  has  the  least  uniform  norm  in  that 
set.  When  the  matrix  A  is  positive  definite  its  eigenvalues  can  be  included  in 
one  single  interval  [a,b]  not  containing  the  origin  and  the  polynomial  having 
least  uniform  norm  in  [a,b]  can  easily  be  expressed  in  terms  of  the  Chebyshev 
polynomials.  This  is  the  essence  of  the  Chebyshev  semi  iterative  method  [5] . 

When  A  is  not  positive  definite,  i.e.  when  it  has  both  positive  and  negative 
eigenvalues,  then  its  speotrum  can  be  included  in  two  disjoint  intervals  [a,b], 
[c,d]  with  b  <  0  <  c  ,  and  the  polynomial  having  having  least  uniform  norm  in 
[a,b]  D  [c.d]  can  no  longer  be  easily  expressed.  Lebedev  [6]  considered  this 
problem  and  gave  an  explicit  solution  to  it  for  the  very  particular  case  when  the 
two  intervals  [a,b]  and  [c.d]  have  the  same  length.  Lebedev  does  not  however  give 
any  suggestion  for  the  general  case  and  de  Boor  and  Rice  [4]  showed  how  to 
prsctically  obtain  the  polynomial  of  minimal  uniform  norm  in  two  disjoint 
intervals.  There  remains  however  several  serious  difficulties  with  the  use  of 
these  minimax  polynomials  which  render  the  method  unattractive  and  unreliable. 
Firstly  the  computation  of  the  minimax  polynomial  can  be  particularly  difficult 
and  time  consuming.  Secondly  once  the  minimax  polynomial  is  available  one  must 
perform  the  Richardson  iteration: 

*i+l  "  xi  +  Ti  ri  (2.5) 

where  the  x  ’»  are  the  inverses  of  the  roots  of  the  minimax  polynomial.  These 
roots  should  therefore  be  computed  before  2.S  is  started  but  this  can  be  done  at  a 
negligible  cost.  The  main  disadvantage  of  2.5,  however,  lies  in  its  unstability. 
This  can  easily  be  understood  by  deriving  from  2.5  the  equivalent  equation  for  the 


residual  vectors 


ri+i“ri  -xl  A  rj  ■  (I  TiA^ri 
which  shows  that 

ri+l  “  (I_Ti+l  A)<1  “  Ti  A  A)  r0  U.6) 

Although  the  final  reaidual  r^  ahould  be  small  in  exact  arithmetic,  the 
intersiediate  residuala  r ^  can  be  very  large  as  is  easily  seen  from  2,6  and  this 
stay  cause  unstability  in  iteration  2.5.  Reordering  of  the  paraaieters  x  night  be 
used  to  achieve  a  better  stability  [1]  although  thia  does  not  seen  to  constitute 
a  definitive  remedy.  Note  also  that  the  computation  of  the  roota  of  the  minimal 
polynomial  can  be  itself  an  unstable  prooesa.  Finally  another  drawback  with  the 
uae  of  the  minimax  polynomial  is  that  if  a  polynomial  of  degree  m  is  used  then  the 
intermediate  vectors  x^  with  i<n  do  not  approximate  the  exact  aolution  x.  Only  the 
last  vector  x^  does  in  theory.  The  implication  of  this  is  that  one  cannot  stop  the 
process  and  teat  for  convergence  before  the  m  steps  are  all  performed.  If  the 
accuracy  is  unsatisfactory  after  the  m  steps,  one  has  the  alternatives  of  either 
recomputinng  another  polynomial  with  different  degree  and  restart  2.5  or  perform 
again  m  atepa  of  iteration  2.5  with  the  same  polynomial. 

It  will  be  seen  that  all  of  the  difficulties  mentioned  above  can  be  avoided 


by  using  polynomials  that  are  small  in  the  two  intervals  in  the  least  squares 
sense  i.e.  polynomials  which  minimise  the  norm  with  respeot  to  some  weight 
function  w(x).  This  brings  up  orthogonal  polynomials  in  two  disjoint  intervals 
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3.  Orthogonal  polynomials  1&  ixft  Dla<olnt  Intervals 

3.1.  Basie  definitions  and  notation 

In  this  subsection  we  will  describe  general  orthogonal  polynomials  in  two 
disjoint  intervals.  Let  (a.b).  (c,d)  be  two  open  intervals  snch  that  b  <  0  < 
c.  The  assumption  that  the  origin  lies  between  the  two  intervals  it  not  essential 
for  this  subsection.  Ve  consider  a  weight  function  w(x)  on  the  interval  (a.d) 
defined  by 


w^(x)  for  x  €  (a.b) 
w(x)  -  w2(x)  for  x  €  (c.d) 

0  for  x  €  (b.e) 

where  w^x)  and  w2(x)  are  two  nonnegative  weight  functions  on  the  intervals  (a.b) 
and  (c.d)  respectively.  The  inner  product  associated  with  w(x)  will  be  denoted  by 
(...).  i.e. : 

<f.g>  “  f (x)g(x) w(x)dx 

or: 

<f»*>  ■  f(x)g(x)w1(x)dx  +  f (x)g(x)w2(x)dx  (3.1) 

The  norm  associated  with  this  inner  product  will  be  denoted  by  II.MC  i.e. 
Ilfllg  -  <f,f>1^2.  Ve  will  often  use  the  notation  <f»g>2  and  <f.g>2  for  the  first 
and  the  second  term  of  the  right  member  of  3.1  respectively. 

The  theory  of  orthogonal  polynomials  shows  that  there  is  a  sequence  (pn)  of 
polynomials  which  is  orthogonal  with  respect  to  <.,.>.  All  of  the  properties  of 
the  orthogonal  polynomials  in  one  interval  hold  for  the  case  of  2  intervals, 
provided  we  consider  these  polynomials  as  being  orthogonal  on  (a.d)  rather  than  on 


(a.b)U(c.d).  For  example  although  ve  oannot  aaaart  that  the  roota  are  all  ia 
[a,b]U[c,d],  we  know  that  they  belong  to  the  interval  [a.dl.  Caaea  where  one  root 
ia  inaide  the  interval  (b.c)  can  easily  be  found.  But  the  fondamental  three  tens 
recurrence  relation  which  holda  for  all  aequencea  of  orthogonal  polynomial a  ia 
obviously  valid.  Thus  the  polynomials  p&  satisfy  a  recurrence  of  the  form: 

Pn+lpn+l(x)  "  (x_an)pn<x)  ~  Pnpn-l(x)  (3*2) 

with 

°n  “  <xpn'pn> 

Pn+1  "  ^pn+l^C  “  <pn+l*pn+l>1/2 

where 

pn+l  <x)  *  (x_an)pn(x)  ~  Pnpn-l<x)  <3*5) 

One  important  problem  considered  in  this  paper  is  to  find  weight  functions 

for  which  the  coefficients  o^,  0^  are  easy  to  generate  and  to  propose  algorithms 

which  generate  them.  Before  entering  into  the  details  of  this  question  in  the 

next  subsection  let  us  briefly  outline  how  these  polynomials  may  be  obtained. 

First  select  two  weight  functions  w^  and  W2  such  that  for  each  of  the  intervals 

(a,b)  and  (c,d)  the  orthogonal  polynomials  considered  separately  are  known.  Then 

express  the  orthogonal  polynomials  in  (a»b)U(c»d)  in  terms  of  the  orthogonal 

polynomials  in  (a,b)  and  in  terms  of  the  orthogonals  polynomil as  in  (c,d).  Ve  thus 

obtain  two  sets  of  n+1  expansion  coefficients  expressing  the  polynomial  p  .  Using 

these  coefficients  we  are  able  to  obtain  the  recurrence  parameters  a  and  0  of 

si  n 

equation  3.2  by  formulas  3.3  and  3.4.  Hence  the  two  sets  of  expansion  coefficients 
can  now  be  updated  by  use  of  3.2  to  yield  the  expansion  coefficients  of  the  next 
polynomial  pa+j  •*••••  Next  the  detaila  of  thia  procedure  will  be  developed  for  a 
judieioua  choice  of  w^  end  w^. 


(3.3) 

(3.4) 
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3.2. 


General lied  Chebyshev  polynomials. 


Let  c1=(a+b) 12,  c^^tc+d) /2  be  the  Biddies  of  the  tvo  intervals  and 
d^fb-a )/2,  d2=(d-c)/2  their  half  widths.  Consider  the  following  weight 
functions: 

w1(x)  -  (2/ it)  Id*-  (x-Cj)2  1"1/2  for  x  €  (a.b)  (3.6) 

w2(x)  *=  (2/n)  [d22-  (x-c2)2  ]~1/2  for  x  6  (c.d)  (3.7) 

The  orthogonal  polynomials  in  (a.b)  with  respect  to  the  weight  function  w^  can  be 
found  by  a  simple  change  of  variables  to  be: 


Tn  1  <3.8) 

where  T&  is  the  Chebyshev  polynomial  of  the  first  kind  of  degree  n.  Likewise  in 
the  second  interval  (c.d)  the  orthogonal  polynomials  with  respect  to  w2  are  the 
polynomials: 


Tn  [<*-c2)/d2  1  (3.9) 

We  are  interested  in  the  polynomials  p  (x)  which  are  orthogonal  in 

n 

(a.b)U(c.d)  with  respect  to  the  weight  function  w(x)  defined  by  w^  in  (a.b)  and  w2 
in  (c.d).  In  order  to  generate  these  polynomials  we  simply  express  them  in  terms 
of  the  polynomials  3.8  in  (a.b)  and  in  terms  of  the  polynomials  3.9  in  (c.d)  as 
follows: 

n  .  . 

P_(x)«  I  r,  T  lU-cJ/d,  ]  (3.10) 

i-0  1  1  11 

V*>«  20Sin)  Til(x_c2>/d2  3  <3'U> 

The  following  proposition  will  enable  oe  to  compute  the  coefficients  un  and 
Pn+1  ’ 


given  the  y's  aad  i*1® 
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P£Ppo»iUpn  3_l1:  Lai  £  fee.  i£i  polynomial  o £  degree  n  having  the  following 
expansions  in  (a.b)  and  (c.d)  : 


n 

p(x)  =  X  y.T  I(x-c1)/d1]  for  x  e  (a,b)  (3.12) 

1=0  1  1  1  1 
n 

p(x)=  l  8.T.[(x-c,)/d,]  for  x  €  (c.d)  (3.13) 

1=0  1 


Set: 


«1  -  2y0  +  2.  r. 

1=1  1 
n-1 

T1  =  2y1y2  +  l  Vw 

1*1 


*2 

T1 


“J  + 


n-l 

2riY2 +  tiyw 


<p.p>  =  o1  +  o2  (3.14) 

<xp.p>  =  Cj  oJt  +  c2  <f2  +  d1x1  +d2t2  (3.15) 

Proof:  Conaider  3.14  firat.  Proa  3.1  and  3.12,  3.13  we  have: 
n  n 

<p.p>=<^YiTi[(x-c1)/d1l,^rlTi[(x-c1)/d1J> 

Using  the  change  of  variable  y  =  (x-c^l/d^  and  the  orthogonality  of  the  Chebysbev 
polynomials  we  obtain 


2  a-  2 

<p  .  p>,  »  2r0  +  l  l .  -  o, 
i=l 

By  the  same  transformations  we  would  obtain  <p,p>2  “  «2.  Bence  <p,p>  =  <p,p>^  + 
<P»P>2  "  aj  +  ®2  Precisely  3.14. 


In  order  to  prove  3.15  consider  <xp,p>2  and  <xp,p>2  separately.  We  have: 


rx~°l] 

<  *P  »  P  >2  -  d2< 


P  >j  ♦  Cj<  p.  p  >2 


(3.16) 


The  first  inner  product  can  be  transformed  as  follows: 
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<  -^p.  P  >x  “  <  1  Tj^ItjU.  2y£  Tji^l]  >! 


,  b  x-c,  x-c,  X-O- 

J  /  2  t4-j-  Ti[_d“  ]*  1  1  dwl< 


x-o, 
1 


x-c, 


where  2  stands  for  2  •  By  the  change  of  variables  y=--r-  this  siaplifies  into 
i-0  ai 


+1 


c^l  p  *  P>!  -  \  /’*  I  I  r4  y  T^y)  )2(l-y2)1/2  dy 


(3.17) 


Using  the  relations: 


y  T^y)  =|  ITi+1(y)  +  T^Cy)  ] 


y  T0(y)  =  T^y) 


we  find 


2  YjT  ^(y)  -  jTj  Tfl(y)  +  <Y0  +  Y2/2)  T2<y> 

+  l(Yi-l  +  W  Ti(y)  +***  +5<Yn-l  +  ^n+l*  Tn(y) 

For  convenience  rn+j  is  set  to  zero.  Replacing  this  in  3.17  and  using  the 

orthogonality  of  the  Chebyshev  polynomials  with  respect  to  the  weight  function 

„  2.1/2 

(1-y  )  we  get 

x_cl  1 

<—f-  P  »  P  >j  =  jp^Yp  *o+*Yl+Y2^2^Yl*l+ 

l(ri+l+Yi-l)ri»i  +***+  l(rn+l  +  Yn-l>Yn*n 
with  s0  -  2  and  Sj^  ■  1  for  i  ^  0  .  This  finally  yields 

<^lp  ,  p  >  -  t1 


and  by  3 .14 


<  *P  .  P  >2  -  OjOj  +  d1r1 


(3.18) 


A  similar  calculation  would  give 
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<*P-P>2  -  c2°2  *  d2X2  (3.39) 

Adding  3.18  and  3.19  finally  yielda  the  deaired  result  3.15. 

Q.E.D. 


Once  the  parameters  and  Pn+i  **e  ooaputed  by  an  appropriate  nse  of  3 .14 

and  3.15,  it  is  possible  to  generate  the  next  orthogonal  polynomial  p  ,,(x).  The 

11+1 

next  proposition  provides  a  simple  way  of  obtaining  Pn+1  when  the  expansions  3.10 
and  3.11  are  used. 


Propoyill.Qfl  3^2:  Let  ^(x) ,  £=0 ,1 , . . .  feg.  the  sequence  of  orthogonal  polynomials 
d&rJLvM  i£SS  1.2.  Sms  lie  expansion  coefficients  y|~)  of  jy^  expansion  3.10 
tAti.sfy  Ihe  relation: 

Pn+lY<iB+1)  =  +  (<sl“0n)^l“>-pnirin"1)'  «.*.•••■+!  <3.20) 

Hill  1M  initial  conditions  : 


B  v(n+1) 
Pn+1Y0 

B  „  <  n+1  ) 

Pn+1T1 


Vi 


*  <vVrJ*)  -  V.™ 


di(Ton)4r2n))  +  -  M 


(n) 

'1  "n'^1 


(n-1) 
n'l 


(3.21) 

(3.22) 


A  re.lg.tl  op  SRPlOSQ.ft*  12  i .JO  gnd  ia  which  il  replaced  fey  *nd  fcy  holds 
1°£  1M  6 'I- 


Proof:  From  3 .2  ve  have 

Wn+l(x)  "  (x~an)  *  Yin>  Til(x“°l)/dl1  -  Pn  I  r|“"1>  ^Kx-c^/dj] 
Using  again  the  change  of  variable  y^x-CjJ/dj  for  convenience  we  get 

*«Vi  ’  ‘V^VS*  l  i'H  V'1  -  K,  i  r'"'1’  T!<t> 

Noticing  that  y  T^y)  -  |(  T^ty)  +  T1+1(y))  when  i  l  1  and  y  T0(y)  -  Tj(y) 


we  obtain: 
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Pn+lpn+l  "  21  ^YV  Ti+l(y>  + 

-K  l  rf°  1  i  <y> 

with  the  notational  convention  T_^»  . 

K*1  »-l  '  dl  «T{m)  *  y2*’«>  Tl<y> 

*  i  I'n-i  *  »[«>  Ti<»>  1  »<«1  -  •„>  i  r‘“’  T,(y> 

-0.  I  rj""1’  1,(7) 

with  the  conventions  stated  in  the  proposition.  Returning  to  the  x  variable  gives 
the  desired  expansion  of  pn+^  In  the  fora  3.10  and  3.11  and  establishes  the 


Vl(y)1  +(cl"an)  i-t  i  Ti(y) 

Therefore 


The  results  of  the  previons  two  propositions  can  now  be  oombined  to  provide 

sn  algorithm  for  coapnting  the  sequences  of  orthogonal  polynomials.  The  ides  of 

the  algorithm  is  quite  simple.  At  a  typical  step  the  polynomial  p^  is  available 

through  its  expansions  3.10  and  3.11  so  one  can  compute  by  formula  3.14.  Then 

we  can  compute  the  expansion  coefficients  of  pn+^  using  proposition  3.2.  Finally 

Pn+1  is  obtained  from  3.4  and  3.14  and  pn+1  is  normalized  by  Pn+1  to  yield  the 

next  polynomial  p  . , .  This  describes  schematically  one  step  of  the  algorithm 
n+1 

given  below.  A  notational  simplification  is  achieved  by  introducing  the  symbol 
2  '  defined  by 

i»n  i“n 

2  '  ».  ■  2  a.  +  2  *4  (3.23) 

i-0  1  0  i-0  1 

ALflflUXBl 

Computation  Ihfi.  generalized  ChCb.YlhgY  PpiThflUitli 
1.  jnlt^allza.  Choose  the  initial  polynomial  Pq(z)  in  such  a  way  that 
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Y0  •"60  1  •  °2  2 


2.  If  rata.  For  n-0,1,2, . . . . ,  a  do 

BAX 

-  Compote  a  by  formala  3. IS: 

11 


(0)  (0)  „ 

T1  :-y2  0 


o  c.o<n)+c  o(n)+d  r(n)  +  d  t(n) 
n  11  22  11  22 


-  Compote  expansion  coefficients  of  p  +.  by  formolas  3.20  -  3.22  and 

tbeir  analogne  for  the  6  ' s .  n 

*  r;:>» ♦  -  v;-» 

i“2  #  3  $  •  •  •  $  ii+l 

sr<a+l)  .  .(m)  .  __  \.<a>  A  .(n-1) 

*0  :“d26l  +(c2"°n)60  "Pn60 

j*(n+l)  /c  (n)  .1  (n)x  . ,  ..(n)  a  .(n-1) 

*1  :=d2(60  )+(c2"an)6l  -pn6l 

r(n+l).3(x^)  +  s<®>)  +  (e  o  ,(n-l) 

*i  (Vl  +  5i+l)  +  (02  °n)Ri  Pn6i 

i*2  »3 1  •  •  •  »ii+l 

-  Compote  llpa+1ll  by  formols  3.14  and  get  0n+1  by  formola  3.4: 


I  1  C?‘“rt>)2 

i-1  1 

J<®+1):«  I1-  [6<»+1>]2 


a  ._  t«(n+1>-te(n+1)  il/2 
pb+i  •  l®i  *°2  J 


-  Normalize  Pji+1  by  0n+1  to  get  p#+1: 


IS 


i  -  0.1,. ...n+1 

-  Compote  new  o's  end  t'i 

(n+1) ,  ~(n+l) ..2  (n+1)  ~(n+l) /ft2 

1  *  1  /Pn+1  °2  2  /Pn+1 

(n+1)  f  ,  (n+1)  (n+1)  (n+l),m  z  ,  ,(n+D  ,(a+l> 

X1  '  £0  Yi  Yi+1  *  X2  *  £0  Ri  6i+l 

At  the  nth  step  of  the  algorithm  one  performs  approximately  20(n+l) 
arithmetic  operations.  Concerning  the  storage  requirements,  we  must  save  the 
coefficients  *•  well  as  the  previous  coefficients  and  Sq11^, 

which  means  4(n  +1)  memory  locations,  if  n  is  the  maximum  degree  allowed. 

This  suggests  that,  in  practice  n^^  should  not  be  too  high.  Note  that  restarting 
could  be  used  if  necessary. 


3.3.  The  minimal  polynomial 


Consider  the  orthonormal  polynomial  p  of  degree  n  produced  by  Algorithm  1. 

According  to  section  2,  an  interesting  approximate  solution  to  the  system  A  x  **  b 

would  be  the  one  for  which  the  residual  vector  r  satisfies: 

SI 

*n  *  Pn<A>  * o 

It  is  easily  verified  that  such  an  approximation  is  given  by 


+  ‘n-l(  A 


)  r. 


where 

s  .(x)  *  [  p  (0)]  1  tp  (0)  -  p  (x)l  /  x  (3.24) 

n-l  u  &  a 

Practically  and  theoretically  however  this  approach  faces  several  drawbacks. 
First,  it  may  happen  that  p  (0)  is  xero  or  very  close  to  zero.  In  that  case  the 

SI 

n**  approximation  x  either  is  not  defined  or  becomes  a  poor  approximation  to  x 
n 
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respectively.  Although  this  difficulty  slight  be  overcome  by  computing  *n+1 
directly  frees  x^,  this  will  not  be  considered. 

A  second  dissdventsge  of  using  the  orthogonel  polynomials  p^  is  that  they  do 
not  sstify  e  suitable  optimality  property.  More  precisely  the  norm  I  lpa/Pn(°) 1 
is  not  minimum  over  all  polynomials  of  degree  not  exceeding  n  satisfying  Pfl(0)  *= 

1.  This  fact  deprives  us  from  any  means  of  comparing  the  residual  norms,  and 
above  all  does  not  ensure  that  the  approximate  solution  x&  will  eventually 
converge  towards  the  exact  solution  x. 

These  remarks  lead  us  to  resorting  to  the  polynomial  %n  of  degree  n  such  that 
q^  (0)  *=  0  which  minimizes  the  norm  llqllg  over  all  polynomials  q  of  degree  £  n. 

Writing  the  polynomials  q  as  q(x)  -  1  -  x  s(x)  whore  s  is  a  polynomial  of 
degree  not  exceeding  n-1,  it  is  clear  that  our  problem  can  be  reformulated  as 
follows: 

e 

Find  a  polynomial  ^  ,  such  that: 

lil-x  s*_j (x)  llc  i  lll-x  s (x) i I c  .  V  s  €  Pa-1  (3.25) 

Here  l-xs(x)  denotes  the  polynomial  q  defined  by  q(x)  »  1-  xp(x)  rather  than 
the  value  of  this  polynomial  at  x. 

From  the  theory  of  approximation  it  is  known  that  the  best  polynomial  s&_^ 

exits  and  is  unique  [3],  Furthermore  the  theory  of  least  squares  approximation 

e 

indicates  that  the  solution  s  ,  must  satisfy  the  equations: 

<l-xs*_1,  x  s(x)  >  -  0,  Vs  €  Pa-1  (3.26) 

To  solve  3.26  let  us  aassume  that  we  can  find  a  sequence  of  polynomials 
q^, j-0,l,..n-l,..  such  that  the  polynomials  xqj(x)  form  an  orthonormal  sequence 
with  respect  to  <  ,  >.  Then  expressing  s*  t  as 


we  obtain  fro*  3.26 


t|j  *  <1.  > 


(3.27) 


It  turns  ont  that  the  polynomials  q^  (x)  can  be  obtained  quite  simply  by  an 
algorithm  which  is  nothing  bat  algorithm  1  with  a  different  initialization.  This 
is  becanse  we  actually  want  to  compute  an  orthogonal  sequence  of  the  form  {xq^} 
with  respect  to  the  same  weiaht  function  w(x)  as  before. 

Practically  we  will  have  again  to  express  the  polynomials  xq^(x)  by  two 
different  formulas  as  in  subsection  3.2: 


n— ^  ' (n) 

xa  (x)  -  It,  t.(x).  x  e  la,b] 
i=0  1  1 

xq  (x)  -  l  6  ,  'n'  T. (x) .  x  e  [c.dl 
^  i-0  1  1 


(3.28) 

(3.29) 


Denoting  by  a*  and  0 *  the  coefficients  of  the  new  three  term  recurrence  which 
n  n 


is  satisfied  by  the  polynomials  q^  we  obtain: 


e;+l  Vl(x)  "  <x_  a;>Vx)  ~Pn^n-l^x^ 


(3.30) 


which  is  equivalent  to: 


Pi+1  x*n+l(x)  "  (x_  x«n<x)  xVl(xl 


(3.31) 


The  condition  that  the  sequence  {  xq.  }  .  .  ,  is  orthonormal  gives  the 

j  j“v  # . .n. . 

following  expressions  for  a'  and  0': 

s  & 


d»<x  (xq^) ,  xq^  >  (3.32) 

Pi+l“  l,x~«n,,C  “  ,,<x  ~  ai)x«n(x)  ’  PixVl(x)  11 C  (3*33) 

It  is  dear  from  the  above  formulas  that  the  Algorithm  for  computing  the 
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sequence  xq  (x)  is  essentially  algorithm  1  in  whioh  the  initialisation  is 
11 

replaced  by  a  new  one  corresponding  to  the  starting  polynomial  xq^  (z)  ■=  z  instead 
of  p^  (z)  -  1.  Thus  the  orthogonal  sequence  xq^(x)  replaces  the  sequence  P&(x)  of 
algorithm  1. 


Ve  must  now  indicate  how  to  obtain  the  expansion  coefficients  q  of  the 

J 

polynomial  s*  in  the  orthogonal  basis  {  q.  } .  According  to  3.27  we  have 
a  4 


qj-<l,zqj> 

Hence 


"j  ■ <  V«Vi  *  <  V‘V» 

-  J(  0.34) 

Here  we  recall  that  the  weight  function  w  has  been  scaled  such  that 

<TQ ,T(j>i«*2,  i«  1,2  and  <T ^ ’ Tj > ^  "  1  for  JdO,  i-1,2.  As  a  consequence  of  3.34  the 

expansion  coefficients  q  can  be  obtained  at  a  negligible  cost.  Ve  can  now  write 

n 

down  the  algorithm  for  computing  the  orthogonal  sequence  qn(x) .  For  simplicity  the 
prime  symbols  will  be  dropped. 

ALGORITHM  £ 

1.  Initialisation.  Choose  the  polynomial  xq^x)  in  such  a  way  that 
qg-constant  and  llxq^llg  ■  1. 


t2:-2(cj  +  c 

2>  +  dj+d2 

t:-t21/2 

ro°):mci/tf 

&0°):“C2/t 

T<0> :“dj/t. 

»{0):-d2/t 

•r- 

+  dj]/t2,  «r*0>:- 

[2c*  +  d2] / t2 , 
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n 


T(0)-  2  Y(0)  V<0)  T(0)-  2  6<0)  6<0) 

T1  *  2  r0  rl  ’  X2  0  ®1  * 


2.  Itmti.  For  n>0,l,2, . . . . ,n  do 


-  Compote  a&  by  formula  3.15: 


v  *  V2*’ 


Compote  expansion  coefficients  of  by  formula  3.20  -  3.22  and 

tbeir  analogne  for  tbe  6  'a. 


-( n+l) 


(n  ,1  (n), 


(n)  a  _ (n-1) 


~(n“?,>  di(  (n)  .  (n).  .  .  (n)  0  (n-1) 

7 1  "21(irj-i  +  ^i+l'  +  'cl“°n>Ti  “  Pnri 


i“2 ,3  .... ,  n+2 


~(n+l)  .(n)  .  x.(a)  o  .U-l) 

80  :  d28l  +<c2"°n)60  Pn60 

~(n+l)  .  /.WA-Wur  ».(n)  a 

81  :md2%  +F2  )+(c2'°n)6l  -Pn6l 

^(n+D.jd-.  (n)  .(nh  +  #.  __  x.(n)  _  R  .(a-D 
6i  :"22(6i-l  +  8i+l)  +  <02  °n)8i  Pn8i 


i-2.3.. . .,n+2 

-  Compote  Uq^jll  *>7  formola  3.14  and  get  0n+1  by  formola  3.4: 

I'"*11:-  l1-  I?’**11]2  .  T •[»<**1>]2 


B  •-  [?<a+1)-ta(n+1>  ]1/2 
Pn+1'  1®1  +a2  1 

Normalise  by  0a+1  to  get 


(n+l)  ~(n+l) 


/ft  »<n+1> .-a<n+l)/ft 

/Pn+1  *  8i  •-8i  /pn+l 


i  ■  0,1 >...,0+1 


-  Compote  new  o' a  and  t’s 
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(n+1)  ~(n+l)  /ft2  (n+1)  ~(n+l)  .  2 

°1  ""l  /Pn+1  °2  '^2  /Pn+1 

(n+1)  f  ,  (n+1)  (n+1)  (n+1) f  ,R(n+l). (n+1) 

i  *  AY*  Y*+1  *  2  *  A  1  1+1 

Note  that  part  2(  iterate  )  ia  identical  with  part  2  of  algorithm  1.  It  is 
iaportant  to  realise  that  the  above  algorithm  does  not  yield  directly  the 
polynoaiials  q  bat  rather  provides  their  expansion  coefficients  in  terms  of  the 
Chebyshev  polynomials  in  the  two  intervals.  In  the  subsequent  sections  ve  will 
need  to  generate  the  polynomial  s  (x) .  This  can  be  achieved  by  coupling  with  part 

A 

2  of  the  above  algorithm  the  following 

Vl(x):"  ri,I<*"tE»)qn<x)  “*n«n<x)1 

A+l 

(n+1)  t(n+l). 

W2<to  +  8o  ) 

•n+l!-  *n  +  VlVl 

starting  with 

qQ(x)-l/t  ,  s0(x)  «  2(yJ0)  +  *J0)  )  qQ(x) 

where  t  is  defined  in  the  initialization  of  the  algorithm. 

3.4.  Compariaon  with  the  minimax  polynomial. 

The  purpose  of  this  subsection  is  to  compare  the  norms  of  the  minimal 
polynomial  described  in  the  previous  subseotion  and  the  polynomial  having  the 
least  uniform  norm  in  [a,b]  U  [o,d]  subject  to  the  constraint  p  (0)  ■  1.  We  will 
denote  by  llpll^  the  uniform  norm  on  [a,b]  U  [c,d],  that  is: 

I  Ip  1 1  -  max  lp(x)l 

x€  [a,b]U(c,d] 

Let  us  denote  by  p*  the  minimal  polynomial  1  -  xs  , (x)  of  degree  n,  defined 

A  A“i 

in  section  3.3,  and  by  p^  the  polynomial  of  degree  n  of  least  uniform  norm  in 


[a,b]U[c,d].  In  other  words  we  have 


(3.35) 


llp*Mc  i  llpllr.  Vp€P, 


I  Ip’ll,  i  llpll..  Vp«p, 


(3.36) 


Since  P  is  a  finite  dimensional  space  all  the  nones  of  P  ere  equivalent  and 
n  n 


the  following  lemma  indicates  precisely  a  comparison  of  the  nones  ll.ll-  and 


the  following  lemma 

ll.ll. 

Lemma  3.3: 

Let  f  1 
n 

1) 

2) 

Proof. 

i) 

n  « 

>1/2 1 


n  C 


(3.37) 

(3.38) 


llfnMi-<f«'  V“<fn'  Vl  +<fn  '  fn>2 

-/  f^(x)  wt(x)  dx  +  /  f^(x)w-(x)  dx 
n  jl  n  z 


i  Ilf  II*  C  /  Wj(x)  dx  +/  w_(x)  dx  ] 

a  o 


and  finally 


ii'X  "'X 


which  gives  the  first  result  by  taking  the  square  root  of  both  members. 
2) 


Ilf  II  -  max  If  (x)l 

n  x6  [a.bJDCc.d]  n 

n  n 

Ilf  II  -max  (max  I  i  r.T.Kx-c.l/d,]  I  ,max  I  I  6  .T.I(x-c,)/d,]  | }  (3 .39) 
*  xCla.bJ  i-0  1  1  1  1  x€[c,d]  i-0  1  1  2  2 
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where  the  y's  end  the  8*s  ere  the  ezpaasioa  ooeff ioients  of  f  In  [a.bl  end  [c.d] 

n 

with  respect  to  the  polynoaials  3.8  end  3.9  respectively.  Using  the  Schwartz 
inequality  we  get 


?  2,1/2 


I  >  T1T.t(x-c1)/d1] I  it  I  rfl 
i=0  1  1  i-0  1 


l  I  T^I(z-c1)/d1]] 
i-0  1  11 


1/2 


Hence 


I  I  y.T  t(*-c  )/d  ]|  i  (n+l)1/2  [  l  y2l1/2 
i-0  1  1  1  1  i-0  1 


Likewise  we  cen  show  that 


,1/2  r  =  .2.1/2 


I  l  6  T  [(z-c  )/d,]|  i  (n+l)i/2  [  l  6‘]] 
i-0  1  1  2  2  i-0  1 

Replacing  3.40  end  3.41  in  3.39  we  get 

Nf  II.  i  (n+l)1/2  >ez{  [  f  y2]1/2.[  I  r ?11/2  } 

i-0  a  i-0  1 

i  <n+l>1/2  [  I  (r 4+6 j)] 1/2 

i-0  1  1 

i  (n+l)1/2  [  I  »  (y2+62)]1/2 
i-0  1  1 

-  (n+l)1/2  ||f  | I 
n  t 


(3.40) 


(3.41) 


Here  we  have  used  the  result  of  proposition  1  end  the  notation  3.24 


Q.E.D 


As  an  immediate  consequence  of  3.37  and  3.38  we  can  easily  bound  each  of  the 
no  rats  II  ll^  and  II  ll^  in  terns  of  the  other  frosi  the  left  and  the  right  as  in 
the  nezt  inequalities: 


(n+l)"1/2  llfjl.i  Illallci2  llfjl.  (3.42) 

\  HV'c  *  ,,fnM-  *  (B+1>1/2  Mfallc  (3.43) 


Ve  now  stats  the  main  result  of  this  subseetioa  which  coapares  with  the  saae 
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•  00 

norm  the  polynomials  p&  and  pQ  . 

CO  • 

Theorem  3 .4:  Let  p  and  p  be  the  minima*  polynomial  and  the  minimal  polynomial 
n  n 

defined  by  3 .35  and  3 .36.  Then  the  followina  inequalities  hold: 

Up*IIc  i  1 1 p » 1 1 c  i  2<tt+1>1/2  llp*llc  (3.44) 

llp"llB(  1  lip*!!.  1  2(n+l)1/2  I  lp“l  1^  (3.45) 

Proof:  Ve  will  only  prove  3.44  as  the  proof  for  3.45  is  identical.  The  first 
part  of  the  double  inequality  is  obyious  by  equation  3.34.  The  second  part  is  a 
simple  consequence  of  lesma  3.3: 

I  IP“I  lc  <  <n+l)1/2  llp“ll„  i(n+l)1/2  NpX 

i  2(n+l) 1/2  llp*llr 
n  t 

Q.E.D. 

•  CD 

The  meaning  of  the  theorem  is  that  the  two  polynomials  p  and  p  have  a 

n  n 

comparable  smallness  in  [a,b]  C  [c,d].  This  will  haye  a  quite  interesting 

implication  for  the  numerical  methods  using  polynomial  iteration  since  it  means 

that  we  can  replace  the  minima*  polynomial  by  the  least  squares  polynomial,  which 

is  by  far  easier  to  generate  and  to  use,  and  still  obtain  a  residual  norm  which  is 

not  larger  than  2(n+l)  times  the  norm  corresponding  to  the  minima*  polynomial. 

Note  that  there  is  no  reason  a  priori  why  to  compare  the  two  polynomials  with  one 

particular  norm  rather  than  the  other.  If  the  norm  II  | |  is  chosen  then  the 
•  00 

polynomial  p  is  smaller  than  the  polynomial  p  while  with  the  uniform  norm  the 
n  n 

contrary  is  true.  Whichever  norm  is  used  however,  the  two  polynomials  have  a  norm 

1/2 

comparable  within  a  factor  not  exceeding  2(n+l)  ,  as  is  asserted  in  the  above 


theorem 
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4.  Application  lit  lftlallfift  al  iftxift  aadtiiaif  aiiam 


4.1.  The  Generalized  Chebyshev  Iteration 


Let  ns  return  to  the  system  2.1  and  suppose  that  ve  have  an  initial  guess 
of  the  solution.  We  seek  for  an  approximate  solution  of  the  form 


x  =  x-  +  s  ,(A)rrt 
n  0  n-1  0 


where  s  ^  is  a  polynomial  of  degree  n-1  and  r^  the  initial  residual.  The 

residual  vector  r  of  x  is  given  by 
n  n 


r  =  [I-As  -(A)]  rft 
o  n-JL  u 


(4.1) 


(4.2) 


The  development  of  the  previous  sections  lead  ns  to  choose  for  sn_j  the  polynomial 

s*  ,  which  minimizes  I |l-xs(x) I over  all  polynomials  s  of  degree  £  n-1.  From 
n-1  t 


subsection  3.3  we  know  that: 


**,(*)  =  I  ^444(*> 
n  1  j-Q  i  i 


(4.3) 


where  the  qj(x)  satisfy  the  recurrence  3.30  and 

■>,  -  *  *J<J)  > 


Let  us  assume  that  the  coefficients  a  ,  8  as  well  as  the  y's  and  6's  are 

n  & 

known.  Noticing  that 

Cl<*)“  •n<x)  +  Vl  Wx) 

and  letting  u^  m  q^(  A  ) rQ ,  it  is  immediately  seen  from  subsection  3.3  that  the 
approximate  solution  xq  can  computed  from  the  recurrence: 


\  •  2<'i 


,u.  . 


x  x  +  q  u 

n+1  n  n  n 
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u  xi o-T  KA-o'  I  )u  -  0'  u  ,  ] 
n+1  0 '  n  n  rn  n-1 


(4.5) 


n+1 


Clearly  the  coefficients  0 '  ,  a'  and  q  can  be  deternined  at  the  sane  tine  as 

n  n  & 

4.5  and  4.4  are  perforned.  Therefore  a  typical  step  of  the  algorithn  would  be  as 
follows: 


ALOOKITOB 

1.  Compote  o'  ,  0 '  .  and  the  coefficients  yt*n+1*,  5(*n+1*  (see  part  2  of 

Algorithm  2).  n  1  1  1 

2.  Compote : 


2[r6U)  +  8o<n)l 


(4.6) 


X  ,  a  *  X  +  II  U 

n+1  n  *n  n 


un+l:=  FLI(A-°;  1  )nn  -  K  an-l  1 


(4.7) 

(4.8) 


n+1 

We  will  refer  to  this  algorithm  as  the  Generalized  Chebvshev  Iteration 

algorithm.  Although  we  need  not  save  the  coefficients  a'  and  0',  if  a  restarting 

n  n 

strategy  is  used  in  which  the  iteration  is  performed  with  the  sane  polynonial  it 
will  be  nore  efficient  to  save  the  coefficients  o'  and  0'. 


Note  that  the  noltiplication  by  a  scalar  in  4.7  can  be  avoided  V.y  using 
instead  of  the  vectors  oq  the  sequence  of  vectors  u'n  defined  by: 

u'  *  n  u 
n  an 

which  can  be  updated  in  no  nore  operations  than  are  needed  for  u  . 

& 

Each  step  of  the  above  algorithm  requires  one  mstrix  by  vector 
multiplication,  3N  additions  3N  multiplications  and  0(n)  operations,  at  the  nt8 


step.  Four  vectors  of  length  N  need  to  be  stored  in  memory  plus  an  extra  0(n) 
locations  for  the  coefficients  y  and  6.  A  feature  of  the  Generalized  Chebyshev 
Iteration,  which  is  particularly  attractive  in  parallel  computation  is  the  absence 
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of  any  inner  product  of  N-dimensional  vectors. 


4*2.  A  minimal  residual  type  implementation 

A  look  at  the  previous  algorithm  reveals  a  similarity  with  one  of  the  various 
versions  of  the  conjugate  gradient  siethod  which  is  described  in  [2] .  One  slight 
ask  whether  it  is  possible  to  provide  a  version  of  the  method  which  resembles  the 
conjugate  gradient  or  the  conjugate  residual  method.  Ve  are  about  to  show  here 
that  the  answer  is  yes  and  we  will  provide  a  method  which  is  very  similar  to  the 
conjugate  residual  algorithm  described  in  [2]. 

To  be  more  specific  we  are  looking  for  an  iteration  of  the  form  : 


*n+l  =  Xn  +  an  wn 

r  .,  ■>  x  -  a  Aw 
n+1  n  n  n 

w  ■  r  . ,  +  b  w 
n+1  n+1  n  n 


(4.9) 

(4.10) 

(4.11) 


where  the  scalars  a  ,  b  are  some  coefficients  that  are  computed  in  such  a  way 
n  n 

that  at  each  step  the  residual  r^+1  is  precisely  the  residual  obtained  by  using 
the  minimal  polynomial  as  described  earlier. 


There  are  two  main  reasons  why  such  a  new  form  of  algorithm  3  might  be  sought 
for.  The  first  is  that  the  above  version  is  slightly  more  economical  than 
algorithm  3,  specifically  4,9  to  4.11  can  be  performed  with  N  less 
multiplications.  Note  also  that  N  less  storage  locations  are  required  in  the  case 
where  the  matrix  by  vector  multiplication  can  be  performed  oomponentwise,  i.e.  if 
the  i-th  component  of  Ax  oan  be  delivered  oheaply  for  all  i.  The  second  reason  for 
looking  for  a  version  of  the  form  given  above  is  that  the  residual  vector  is 
available  at  each  step  of  the  iteration. 


A 


Some  theoritical  background  is  necessary  for  the  understanding  of  the 

remainder  of  this  subsection.  The  development  of  the  conjugate  gradient  like 

aiethod  to  be  presented  is  based  on  the  simple  remark  that  the  approximate  solution 

x^  belongs  to  the  affine  subspace  Xq+Ko  where  is  the  so  called  Krylov  subspace 

2  H“"l 

spanned  by  r^,  Ar^ ,  ArQ , . . .  A  rQ  .  Then  the  key  observation  is  that  there  is 

a  one  to  one  mapping  between  the  subspace  K  and  the  space  P  .  of  polynomials  of 

n  n-i 

degree  not  exceeding  n-1.  In  effect  for  any  element  u  of  K&  of  the  form  u  -  (ijr^ 

+  +  ....+  lf0  *  we  can  **,ociate  the  polynomial  q^(t)  •  +  (ijt  + 

. ...+|in_ltn  ^ .  Clearly  the  transformation  which  maps  u  into  q^  is  an  isomorphism 

between  K  and  P  , .  Ve  now  introduce  a  new  inner  produot  on  K  . 

n  n-1  n 

Definition:  Ve  will  call  Generalized  Chebyshev  inner  product  on  K^  and  will  denote 
by  (.,.)c  the  bilinear  form  defined  by  : 

(u,v)c  =>  <qB#qy>  (4.12) 

where  <.,.>  is  the  the  inner  product  ob  P  ,  defined  by  3.1,  3.6  and  3.7. 

n-i 

That  this  constitutes  an  inner  product  is  an  immediate  consequence  of  the 
isomorphism  introduced  above.  There  is  no  ambiguity  in  denoting  again  by  ll.ll^, 
the  norm  on  K^  induced  by  this  inner  product  .  Ve  then  have  the  following 
immediate  result. 

Pr.9P9.»it*9n  UL:  -M  ®»ch  jjjp  ai  ite  gmaiiiafl  lianOas*  iM 

aamalaiia  lalatiaa  \  1M  SbaluMJuq  a ioou  I  If-Ax  \\csl 

lb£  St JiflJUi  Ittl fill  f-Ax  over  ill  vectors  2  si  Kn- 

In  other  words  we  have: 

Ilf  -  A  x  ||  i  ||f  -  Ax  ||-  V  x  €  K 
n  t  n 


The  proof  of  the  proposition  is  straightforward 


w 


» 
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The  most  convenient  way  of  deriving  a  veraion  of  algorithm  3  equivalent  to 
4.9-4.11  ia  to  nse  a  little  analogy.  Indeed  ainee  we  want  to  minimize  the  residual 
norm  with  respect  to  the  generalized  Chebyshev  inner  product  4.12  we  can  uae  an 
adapted  version  of  the  minimal  residual  method  which  achieves  the  same  goal  with  a 
different  inner  product.  The  version  given  below  is  precisely  an  adaptation  of  the 
conjugate  residual  method  described  in  [2] . 


jMflBIl  £ 

1.  Initialize.  Choose  an  initial  vector  xfi  and  compute  rQ“f-A*0,  w0=r0< 

2.  ItPXLt*. 


*n  <x*r  *qr  >/<xSr  > 

n  n  n  n 

(4.13) 

x  ..  :*  x  +  a  w 
n+1  n  n  n 

(4.14) 

r  :*  r  -  a  Aw 

n+1  n  n  n 

(4.15) 

b  <xq  ,q  >/<xq  ,a  > 

n+1  n+1  n  n 

(4.16) 

Wn+1  ‘  rn+l  +  **nwn 

(4.17) 

Sr  «r  +  b»S 

n+1  n+1  n 

(4.18) 

The  computation  of  the  inner  products  involved  in  4.13  and  4.16  can  be 

carried  out  in  a  similar  way  as  before  by  use  of  proposition  3.1  together  with 

expressions  similar  to  3.12,  3.13  for  each  of  the  polynomials  q  and  q^  but  we 

rn  n 

omit  the  details. 


As  with  the  conjugate  residual  method  in  the  indefinite  case,  the  above 

algorithm  faces  s  risk  of  breakdown  in  equation  4.16  because  <xq  ,q  >  can  be 

rn  rn 

equal  to  zero.  We  do  not  know  under  what  circumstances  the  inner  product 

<xq  ,q  >  vanishes  but  our  short  experience  is  that  this  does  not  often  happen, 
r  r 
n  n 


We  will  prefer  algorithm  3  in  general  de.pite  the  little  extra  work  and  possibly 
storage  involved.  Note  that  it  is  also  possible  to  write  a  conjugate  gradient  type 
method  which  would  correspond  to  using  as  residual  polynomial  the  orthogonal 
polynomial  of  subsection  3.2,  instead  of  the  minimal  polynomial  of  subsection  3.3. 

4.3.  A  Block  Generalization 

A  block  algorithm  can  be  derived  from  the  Generalized  Chebyshev  Iteration 
described  above.  The  principle  of  the  Blook  algorithm  is  similar  to  that  of  the 
Block  Stieffel  iteration  proposed  in  [11] . 

As  shown  in  fig.  4-1  if  an  eigenvalue  lies  inside  the  interval  [b,c]  then 

e  e 

p  (X.)  is  large  in  comparison  with  the  other  p  (X^),  with  j  ^  i. 

As  a  consequence  the  residual  r^  »  pa(A)rQ  will  be  mainly  in  the  direction  of 

the  eigenvector  z^  asociated  with  X^  Therefore  by  a  projection  process  the 

component  of  the  residual  in  this  direction  can  be  eliminated  leaving  only  the 

small  components  in  the  directions  Zj  with  j  /  i.  Let  us  now  assume  that  instead 

of  one  eigenvalue  X^,  we  had  m  eigenvalues  X^,X^|^ , .... ,X^)m  ^  inside  [b,c],  and 

that  we  are  performing  m  different  Generalized  Chebyshev  iteration,  using  m 

different  starting  vectors  x_  ,xn  ...  .  ,x.  .  We  can  again  make  the  same 

argument:  the  residuals  r  .  at  the  n**1  step  will  have  their  dominant  componants 

n,  x 

in  the  eigenvectors  z  zi+j , . . . ,  and  we  oan  eliminate  these  components  by 

projecting  onto  the  subspace  spanned  by  Cr  ,} ,  ,  ,  .  It  was  shown  in  [11] 

n,j 

that  in  the  positive  definite  esse,  the  rate  of  convergence  of  the  Block  Stiefel 
Iteration  is  in  general  substantially  superior  to  that  of  the  single  vector 
Stiefel  Iteration. 

The  block  algorithm  is  more  particularly  suitable  for  parallel  computation 


since  the  m  Generalized  Chebyshev  iterations  involved  s,  can  be  performed 
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simultaneously.  Ignoring  the  matrix  by  vector  multiplication,  intercommunication 
between  the  processora  ia  needed  only  when  the  projection  proceaa  takes  place, 
which  is  infrequent.  As  mentioned  earlier  the  fact  that  there  are  no  inner 
products  of  N  dimensional  vectors  is  a  further  advantage  of  this  algorithm  in 
parallel  computation. 

4.4.  A  Priori  Error  Bowmds. 

We  now  establish  an  a  priori  error  bound  which  will  show  that  the  Generalized 
Chebyshev  Iteration  is  convergent.  We  should  emphasize  however  that  the  bound 
given  here  is  proposed  for  the  sole  purpose  of  demonstrating  the  global 
convergence  of  the  process  and  does  not  constitute  a  aharp  error  bound  in  general. 

The  next  theorem  is  based  on  the  following  estimate  of  the  residual  norm 
which  is  an  immediate  consequnece  of  lemma  3.3 

Lemma  4.2:  At  the  n-th  step  of  of  the  Generalized  Chebyshev  iteration  the 

residual  norm  r  of  the  approximate  solution  x  satisfies: 
ii  n 

Mr  II  <  (n+l)1/2  II  p*llr  II  rJI  (4.19) 

n  n  t  u 

Proof :  From  equation  2.4  we  have 

Mr  II  1  I  lr-1 1  max  Ip  (A  )  f 
“  °  i«l,..N  n  1 

Therefore 

i  "r0M  "01- 

The  result  follows  from  inequality  3.37  of  lemma  3.3. 

Q.E.D. 


Xhfiju&a  £jl:  Saamm  that  ife  spectrum  a!  A  ii  gsat.ilflid  ia  XjufclBLufli  lad  lai 
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m  *  mint  lb  I  »  I c I ) ,  M  -  msxMal.ldl) 


Then  at  .the  nth  ilfig  q£  the  GpnprjJLjfcfttf  Chebvshcv  Iteration  1M  residual  Tec tor  r 

n 

satisfies  the  ineaoelitv: 


MrJI  i 


2(n+l) 


1/2 


llr0ll 


Tn*[ 


M2-.2 


(4.20) 


where  n*  is  the  integer  division  of  n  tv  2  end  represents  the  Chebvshev 
polynomial  of  i  of  1^2.  first  fcfftfli 

Proof :  By  equation  4.19  we  have 

llrjl  i  (n+l)1/2  ||p*Mc  ||r0ll  (4.21) 

i  (n+l)1/2  llpllc  l(r0lf 

for  any  polynoaial  p  of  degree  f  n  such  that  p(0)  =  1.  Consider  the  polynoaial: 

t  (*)  -  T  ,lq(*)l  /  T  , t q(0) ] 
n  n  n 

with 

q(g)  -  [M2  +  a2  -  x2  ] /IM2  -  b21 

The  degree  of  t  is  <  n  and  we  have  t  (0)  «=  1.  Therefore  froa  4.21  end  leans  3.3 
n  & 

we  get 

llrjl  1  (n+l)1/2  l|tnllc  ||r0ll  i  2<n+l)1/2  I  ItJI  JI*0I  I 

To  ooaplete  the  proof  we  observe  that  t  has  been  chosen  in  snch  a  way  that 

SI 

lit  ||  -[T  ( q(0) )  ]  2  because  for  x  belonging  to  [a,b]U[e»d] »  It  (x)|  £  1,  the 

ii  ®  n  n 

value  1  being  reached  for  sone  x's  in  [a.blOIc.d]  . 


Q.E.D. 


Note  that  ideally  M2  and  m2  are  the  largest  and  the  smallest  eigenvalues  of 

2 

A  .  Hence  an  interpretation  of  4.20  is  that  n  steps  of  our  algorithm  will  yield  a 

1/2 

residual  nora  not  larger  than  2(n+l)  tiaea  the  residual  nora  obtained  froa 

fn/2f  steps  of  the  classical  chebyshev  iteration  applied  to  the  noraal  equations 
2 

A  *  =  A  b.  This  result  is  however  a  pessimistic  estimate  as  the  numerical 
experiments  show. 

5.  Application  £o  computation  al  ciaenwectora  associated  with  interior 
eiaenvalues. 

Consider  the  eigenvalue  problem 

A  x  -  A  z  (5.1) 

where  A  is  N  dimensional  and  symmetric.  When  the  desired  eigenvalue  X.  is  one  of 
the  few  extreme  eigenvalues  of  A,  there  are  a  host  of  methods  for  solving  5.1, 
among  which  the  Lanczos  algorithm  appears  as  one  of  the  most  powerful.  When  X  is 
in  the  interior  of  the  spectrum  then  problem  5.1  becomes  more  difficult  to  solve 
numerically.  The  Lanczos  algorithm  often  suffers  from  a  slow  rate  of  convergence 
in  that  case  and  may  require  a  number  of  steps  as  high  as  several  times  the 
dimension  N  of  A.  A  further  disadvantage  of  this  is  that  the  Lanczos  vectors  must 
be  kept  in  secondary  storage  if  the  eigenvectors  are  wanted.  This  might  be 
acceptable  if  a  large  number  of  eigenvalues  are  needed  but  when  one  only  seeks  for 
a  small  number  of  them,  possibly  with  a  moderate  accuracy,  then  the  use  of 
techniques  based  upon  subspace  iteration  should  not  be  discarded. 

We  would  like  to  indicate  in  this  seotion  how  the  polynomials  described  in 
section  3  can  effectively  be  used  to  obtain  approximations  to  eigenvectors 
associated  with  interior  eigenvalues.  Consider  an  approximation  to  the  eigenvector 
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vn  =  pn  (A)  v0 


(5.2) 


where  v_  is  some  initial  approximation  of  z  and  p  a  polynomial  of  degree 
u  & 

n.  Writing  v.  in  the  eigenbasis  {z.}.  .  „  at 

u  i  i*i 9  * • w 


*0  -  i«i*i 

1=1 


we  obtain: 


vn  ’  i«i  »n(V  zi 
i=l 


(5.3) 


which  shows  that  if  vq  is  to  be  a  good  approximation  of  the  eigenvector  then 

every  p  (X.)  with  j  t  i  must  be  small  in  comparison  with  p  (X.)  .  This  leads  to 
n  j  n  l 

seek  for  a  polynomial  which  is  small  in  '^i-l^^i-t-l  an<*  which  satisfies 

Pa(X^)  =1.  In  other  words  we  need  to  find  a  polynomial  of  the  form 

l-(x-X,)t  , ( x-X . )  which  is  small  in  these  intervals  in  the  least  squares  sense, 

i  n-i  i 

The  simple  change  of  variables  t  *  x-X^  transforms  this  problem  into  the  one  of 
subsection  3.3  thus  providing  an  algorithm  for  computing  an  approximation  to  the 
desired  eigenvector  associated  with  Xj.  Note  that  in  practice  X^  is  not  available 
and  is  replaced  by  some  approximation  p  which  is  improved  during  the  process.  We 
will  find  a  polynomial  in  the  variable  t  which  corresponds  to  the  operator  A  -  pi 
i.e.  the  approximate  eigenvector  will  have  the  form 

Vn  'B  pn  v0  “  <  1  "  ^W1)  ‘^(A-pI)  J  vQ 

where  p^is  the  least  squares  polynomial  obtained  from  algorithm  2  with  : 

a  =  Xx  -  p,  b=  X--p,  c  =  X+  -  p,  d=  XN-p 


where  X^^  and  X^+^  are  replaced  by  X  and  X  respectively. 


An  interesting  problem  which  arises  here  is  to  estimate  the  interior 
parameters  b  and  c,  or  more  precisely  to  refine  dynamically  a  given  pair  of 


starting  parameters  b,c.  This  is  discussed  in  the  next  section  and  a  few 
suggestions  are  given  there. 
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An  obvious  extension  of  this  algortiha  which  we  omit  to  consider  here  is  the 
subspace  iteration  otherwise  known  as  simultaneous  iteration  method  (  See 
e.g.  [9]) in  which  one  iterates  with  a  block  of  vectors  simultaneously  and  uses  the 
Rayleigh  Ritz  projection  approximations  from  the  subspace  spanned  by  the  blocks. 


6.  Some  pmt  iRll  aJUUAiaXMliOUL 


6.1.  Projection  techniques  and  the  estimation  of  the  optimal  parameters. 

Little  was  said  in  the  previous  sections  about  the  determination  of  the 
parameters  a,b.c,d.  It  is  not  often  the  case  that  these  parameters  are  known 
beforehand,  and  one  has  to  provide  means  for  computing  them  dynamically.  This  part 
of  the  algorithm  is  a  determinant  factors  of  the  efficiency.  Let  us  recall  that 
a,d  are  ideally  the  smallest  and  largest  eigenvalues  Xj  and  X^  of  A,  while  b  and  c 
should  be  equal  to  the  eigenvalues  X  and  X+  closest  to  the  origin  with  X  <0  and 
X+>0.  The  smallest  and  largest  eigenvalues  X^  and  X^  are  easier  to  estimate  and 
we  have  several  ways  of  doing  so,  the  simplest  being  perhaps  the  use  of  Gershgorin 
disks.  The  use  of  the  Lanczos  algorithm  is  also  particularly  interesting  since  it 
yields  at  the  same  time  a  good  estimate  of  the  parameters  a,d  and  a  fairly 
accurate  starting  vector  for  the  iteration(see  e.g.  [7]).  It  is  a  fact  that  the 
extremal  parameters  a.d  which  are  easier  to  oompute  are  also  more  important  for 
the  convergence:  if  a  is  larger  than  X^  or  d  is  less  than 
risk  that  the  process  will  diverge.  This  is  because  outside  the  interval  Ia,d)  the 
residual  polynomial  can  have  huge  values  as  is  seen  in  figure  6-1.  When  this 
happens  however  nothing  is  lost  because  we  can  stop  and  use  the  (huge)  residuals 
as  approximate  eigenvectors  associated  with  X^  and  X^  thus  providing  more  accurate 


Xjj  then  there  is  a  high 
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estimates  of  both  a  and  d,  with  the  help  of  the  Rayleigh  quotient.  (See  figure 
6-1) 

For  the  interior  parameters  b  and  c  a  projection  technique  based  on  a  similar 
principle  can  be  developed.  In  effect  if  we  start  with  two  approximations  b  and  c 
such  that  b  <  X  and  c  >  X+,then  as  is  seen  from  figure  4-1  the  residual  vectors 
will  have  large  components  in  the  eigenvectors  associated  with  the  eigenvalues  X 
and  X+  closest  to  the  origin.  A  Rayleigh  Ritz  approximation  can  therefore  by 
applied  to  extract  more  accurate  estimates  of  X  and  X*  from  two  or  more  residual 
vectors. 

In  fact  a  more  complete  procedure  can  be  described  as  follows.  Suppose  that 

we  start  with  an  interval  [b.c]  sufficiently  wide  to  contain  the  eigenvalues  X+, 

and  X  .  Then  after  a  number  of  steps  the  residual  vectors  will  give  an  accurate 

representation  of  the  eigenspace  corresponding  to  the  eigenvalues  contained  in  the 

interval  [X  ,X  ]  .  Therefore  a  projection  process  onto  the  subspace  spanned  by  a 

small  number  of  successive  residual  vectors  will  provide  not  only  a  good 

approximation  of  the  eigenvalues  X  and  X+,  but  also  the  Rayleigh  Ritz  approximate 

solution  will  be  more  accurate  than  the  current  approximation.  In  practice  we 

will  iterate  with  a  polynomial  of  fixed  degree.  Ve  will  orthonormal ize  the  m 

latest  residuals  where  m  is  usually  5  or  6  (in  order  to  have  an  accurate 

representation  of  at  least  three  eigenvectors  ).  Calling  Q  the  orthonormal  set  of 

T 

vectors  obtained  from  this  orthogonal ization,  we  compute  the  Ritz  matrix  Q  A  Q 
and  its  eigenvalues  by  e.g.  the  OR  algorithm.  Ve  then  compare  tA»  Ritz 
eigenvalues  thus  obtained  with  those  of  the  previous  projection.  After  a  certain 
number  of  iterations  some  of  the  Ritz  eigenvalues  will  converge  to  a  moderate 
accuracy  (say  to  2  or  3  digits  of  accuracy).  As  soon  as  an  eigenvalue  which 
corresponds  to  either  X~,  or  X*  converges,  we  replace  b  or  c  by  that  eigenvalue. 
When  both  of  the  eigenvalues  have  converged,  we  replace  the  approximate  solution 
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__  ^  y 

by  the  Rayleigh  Ritz  approximation  x  =  x  +  £  IX.l  z,r  where  the  sum  is  over  the 

n  j  j  n 

converged  Ritz  eigenvalues  A.  and  the  z.  are  the  associated  approximate  Ritz 

J  J 

eigenvectors.  Recall  that  z^  is  defined  by  z^  =  Q.s^  where  s^  is  an  eigenvector  of 
X 

Q  A  Q  associated  with  the  Ritz  value  A^ .  The  effect  of  this  important  projection 
step  is  to  eliminate  the  (large)  components  of  the  residual  in  the  direction  of 


the  converged  eigenvectors  and  is  quite  effective  as  will  be  seen  in  some 
experiments  described  in  section  7.  A  useful  variation  here  is  to  use  the  latest 


directions  of  search  u^  instead  of  the  latest  residuals.  We  thus  avoid  to  compute 

residual  vectors.  The  reason  for  this  is  that  q  (A)r_  is  a  vector  having  large 

n  u 

components  in  the  directions  of  the  eigenvectors  associated  with  the  eigenvalues 

nearest  to  zero.  Our  experience  is  that  the  results  provided  by  the  use  of  u^ 

instead  of  r  are  often  slightly  better  than  those  using  the  uneconomical 
n 

residuals. 


The  above  process  is  even  more  efficiently  implemented  in  a  Block  version, 
because  then  we  have  a  good  approximation  to  several  eigenvalues  at  the  same  time. 
Note  that  in  this  case  we  need  to  have  the  interval  [b,c]  enclose  at  least  m-1 
eigenvalues  if  m  is  the  dimension  of  the  blocks. 

The  same  projection  technique  as  the  one  described  above  can  be  implemented 
for  the  problem  of  computing  an  interior  eigenvalue  A^  and  its  corresponding 
eigenvector.  We  will  drop  the  subscript  i  in  the  following  discussion.  We  are 
again  assuming  here  that  we  already  know  a  good  estimate  of  the  extreme  parameters 
a,  d  and  that  we  start  with  an  interval  Ib,c]  sufficiently  wide  to  contain  the 
eigenvalues  A,  A+,  and  A  .  Then  after  a  number  of  steps  of  the  iteration 
described  in  section  5  the  approximate  eigenvectors  will  give  an  accurate 
representation  of  the  eigenspace  corresponding  to  the  eigenvalues  contained  in  the 
interval  [b,c].  Hence  a  projection  process  onto  the  subspsce  spanned  by  a  small 
number  of  successive  approximate  eigenvectors  will  provide  good  approximations  to 
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the  eigenvalues  X,  X  and  X*.  Also  the  Rita  eigenvector  will  be  a  much  better 

approximation  than  the  current  approximation.  Practically  we  will  orthonormal ize 

the  m  latest  approximations  into  the  N  x  m  orthonormal  system  Q.  Ve  then  compote 
T 

the  Ritz  matrix  Q  A  Q  and  its  eigenvalues  by  e.g.  the  QR  algorithm  and  compare 

the  Ritz  eigenvalues  obtained  with  those  of  the  previous  projection.  As  soon  as 

an  eigenvalue  which  corresponds  to  either  X,  X  ,  or  X+  starts  converging,  we 

replace  |i>  b  or  c  by  that  eigenvalue.  When  the  three  of  the  eigenvalues  have 

converged,  we  replace  the  approximate  eigenvector  by  the  Ritz  vector  z  =  Q  s, 

T 

where  s  is  the  eigenvector  of  Q  A  Q  associated  with  X.  This  process  is  tested  in 
an  experiment  in  section  7. 

6.2.  Using  high  degree  polynomials 

One  might  ask  whether  it  is  possible  to  use  high  degree  polynomials  in 

practice.  Our  experience  is  that  despite  the  fact  that  we  often  encounter 

underflow  situations  with  high  degree  polynomials,  if  these  underflows  are  handled 

simply  by  replacing  the  underflowing  numbers  by  zeroes,  it  is  always  better  to  use 

a  high  degree  polynomial  than  a  low  degree  polynomial.  This  fact  will  be 

illustrated  in  an  experiment  in  the  next  section.  Polynomials  of  degree  as  high 

as  200  or  300  are  quite  useful  for  badly  conditioned  problems.  We  open  a 

parenthesis  here  to  point  out  the  following  interesting  observation.  Jt  has  been 

observed  during  the  numerical  experiments  that  the  last  coefficients  6^, 

become  tiny  as  1  and  n  increase.  This  is  the  cause  of  the  underflowing  conditions 

mentioned  above.  The  practical  consequences  of  this  observation  is  that  after  a 

certain  number  of  steps  we  do  not  need  to  save  those  tiny  elements.  For  high 

degree  polynomials  this  will  result  in  a  nonnegl igible  cut  off  in  memory  needs. 

In  fact  a  more  important  observation  is  that  the  coefficients  a  and  p  are 

n  n 

cyclically  converging  more  precisely  it  seems  that  1*  •  converging  sequence 

for  fixed  j.  The  same  phenomenon  is  true  for  the  p's.  A  proof  of  this  phenomenon 
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is  not  available  however.  As  a  consequence  there  is  a  hope  that  after  a  certain 
step  we  can  simply  use  the  previous  o's  and  p's  thus  avoiding  farther  calculations 
of  these  cofficients. 


7.  Numerical  experiments 


In  this  section  we  will  describe  a  few  numerical  experiments  and  give  some 

more  details  on  the  pactical  use  of  the  Generalized  Chebyshev  Iteration.  All  the 

experiments  have  been  performed  on  a  DEC  20,  nsing  donble  precision  (nnit  roundoff 
-17 

of  order  10  )  .  Any  mention  to  the  Generalized  Chebyshev  iteration  refers  to 

algortihm  3  rather  than  algorithm  4. 

7.1.  Coaqparisom  with  8TMMLQ 

The  SYMMLQ  algorithm  described  in  [8]  and  the  various  versions  of  it  such  as 
MINRES  [8]  SYMMBK  [2]  all  based  on  the  Lanczos  algorithm,  are  very  probably  the 
best  known  iterative  methods  for  solving  large  indefinite  systems  at  the  present 
time  so  it  is  important  to  compare  the  performances  of  our  algorithm  with  this 
class  of  methods.  We  will  mainly  consider  SYMMLQ  although  some  of  the  other 
versions  are  slightly  faster  (but  also  slightly  less  stable).  Although  it  is 
difficult  to  make  any  definitive  conclusion  we  will  see  that  there  are  instances 
where  the  Generalized  Chebyshev  iteration  has  a  better  performance  than  SYMMLQ. 
Table  1  shows  the  work  per  iteration  and  the  storage  requirements  of  both 
algorithms  when  the  number  of  nonzero  elements  in  the  matrix  A  is  equal  to  NZ  (We 
have  not  counted  the  storage  required  for  the  matrix  A) .  An  obvious  advantage  of 
the  Generalized  Chebyshev  iteration  is  its  low  storage  requirement.  Note  that  with 
algorithm  4  GCI  would  require  2N+-NZ  operations  instead  of  3 V+NZ. 


H - 4 

1  Add/Mult- s  I 

Storage 

-+ 

1 

1  SYMMLQ 

I  9  N  +  NZ  | 

6  N 

1 

1  GCI 

j - 

1  3  N  +  NZ  I 

-i - - •» 

4  N 

1 

-  + 

Let  us  discuss  the  above  table  under  the  assumption  that  NZ  ■  5  N,  which 


would  correspond  for  example  to  applying  the  Inverse  iteration  method  to  compute 
an  eigenvector  of  a  Laplace  operator.  Then  it  is  easily  seen  that  our  algorithm 
becomes  competitive  if  the  number  steps  for  GCI  does  not  exceed  1.75  times  the 
number  of  steps  required  by  SYXMLQ.  It  is  observed  in  practice  that  this  ratio  is 
seldom  reached.  In  fact  as  the  next  experiment  will  show,  when  the  problem  is  well 
conditioned  i.e.  when  both  b  and  c  are  not  too  small  relative  to  a  and  d,  then  the 
number  of  steps  is  not  too  different  for  the  two  algorithms.  For  the  first 
experiment  we  have  chosen  a  diagonal  matrix  diagd^)  where  the  eigenvalues  are 
distributed  uniformely  in  each  interval  [a,b]  and  [c,dl.  We  have  taken  N*20t, 
a— 2.0,  b=-0.5,  c**0.5  and  d*=6.0.  The  eigenvalues  X^  are  defined  by: 

define  ij  *  T  N(b-a)/(b-a  +  d-c) 1 

for  i»l ,2 . . ij :  Xj,  -  a  +  (i-D.hj,  with  hj-tb-a)^  ij-l)  (7.1) 

for  i»i1+l....N:  X£  -  c  ♦  (i-D.hj,  with  hj-U-c) /(N-ij)  (7.2) 

T 

The  right  hand  side  f  has  been  taken  equal  to  f  ■=  A  e,  where  e  *(1,1 ,1 . . .1) . 
The  initial  vector  was  a  random  vector,  the  same  for  both  SYMMLQ  and  GCI.  The 
Generalised  Chebyshev  iteration  using  the  exaot  parameters  a,b,c,d  is  run  with  a 
residual  polynoomial  of  degree  25.  The  residual  norms  are  plotted  in  figure  7-1 
in  logarithmic  scale. This  is  done  every  5  steps  for  a  total  of  75  steps.  Observe 
that  the  behaviors  of  both  algorithms  are  almost  identical. 

This  example  shows  the  following  interesting  property:  when  the  origin  is 
well  separated  from  the  spectrum  the  C-G  like  methods  will  converge  with  about  the 
same  speed,  but  each  step  is  more  expensive  as  shown  in  table  1. 

The  next  test  compares  the  two  algorithms  in  presence  of  a  mildly  badly 
conditioned  problem.  Ve  have  taken  the  same  example  as  above  but  the  values  of  b 
and  c  have  been  changed  into  -0.05  and  0.05  respectively,  snd  N  has  been  increased 


44 


to  500.  The  eigenvalues  X^  ere  defined  again  by  7.1  and  7.2.  Tbit  it  sore 
difficult  to  tolve  than  our  previout  example.  Ve  have  plotted  in  fig. 7-2  the 
logarithms  of  the  retidual  norma  for  a  total  number  of  etept  of  500.  The  inner 
loop  of  6CI  utes  a  polynomial  of  degree  250.  (therefore  2  outer  iteration!  were 
taken.)  The  initial  vector  hat  been  choten  randomly  and  hat  an  initial  retidual 
norm  of. 719  E+02.  Although  SYMMLQ  requiret  lett  steps  here  to  achieve  a  reduction 
of  the  residual  norm  by  an  order  of  l.E-06,  observe  that  in  the  beginning  the  GCI 
has  a  better  performance  than  SYMMLQ.  But  at  the  number  of  steps  approaches  the 
dimension  N  of  A,  we  observe  a  speeding  up  of  SYMMLQ.  Note  that  in  exact 
arithmectic  SYMMLQ  would  have  converged  exactly  at  the  500-th  steps.  As  a 
comparison  we  plotted  on  the  same  figure  the  residual  norms  obtained  with  GCI 
using  a  low  degree  polynomial.  Observe  the  slowing  down  of  GCI  after  130  steps  due 
to  the  fact  that  the  residuals  start  being  in  the  direction  of  some  eigenvectors. 

This  shows  the  superiority  of  using  a  higher  degree  polynomial  when  possible. 

In  the  above  experiments  we  have  assumed  that  the  optimal  parameters  a,b,c,d 
are  available.  This  is  unrealistic  in  practice  unless  another  linear  system  has 
already  been  solved  with  the  same  matrix  A.  Ve  would  like  next  to  show  the 
effectiveness  of  the  projection  procedure  described  in  section  6.  In  the  same 
example  of  dimension  500  as  above  we  have  taken  as  extremal  parameters  the  exact 
values  a  and  d.  The  interior  parameters  b  and  c  were  initially  set  to  -0.15  and 
0.15  respectively,  instead  of  -0.05  and  0.05.  Ve  iterated  with  a  polynomial  of 
degree  50.  After  each  (outer  )  iteration  a  Rayleigh  Rita  step  using  the  10  latest 
vectors  n^  as  suggested  in  section  6  was  taken.  The  convergence  of  the  Rits 
eigenvalues  is  very  slow  in  this  example  because  the  relative  gap  between  the 
eigenvalues  of  A  is  quite  small.  At  the  9-th  iteration  no  Ritx  value  has 

_3 

converged  to  the  demanded  accuraoy  10  and  it  was  decided  to  take  a  projection 
step  anyway  with  those  eigenelements  which  have  converged  to  a  relative  accuracy 
of  0.2.  Tbs  code  then  gave  as  eigenvalues  X-— 0.0504...  and  X+»0.0507  and  the 
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final  (10-th)  iteration  was  pursued  with  these  values  for  b  and  c.  Figure  7-3  is  a 
plot  of  the  residual  noras  obtained  with  this  technique  together  with  those 
previously  obtained  by  using  the  exact  paraaeters  with  a  residual  polynoaial  of 
degree  2S0.  It  is  reaarksble  that  the  final  results  are  not  too  far  from  those 
using  the  optiaal  paraaeters. 

7.2.  Coaputing  interior  eigeneleaenta. 

Consider  the  following  N  x  N  five-point  discretization  aatrix  : 
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By  Gershgorin's  theorem  all  the  eigenvalues  of  A  lie  in  the  interval  [0.8]. 

Suppose  that  we  would  like  to  coapute  the  eigenvalue  nearest  to  the  value  2.5  and 

its  corresponding  eigenvector.  Assuming  that  we  do  not  have  at  our  disposal  any 

estimate  of  the  eigenvalues  nearest  to  2.5.  we  start  the  iteration  described  in 

section  5  with  p  equal  precisely  to  2.5.  Also  as  the  first  estiaates  of  the 

eigenvalues  X*  and  X  .  the  saallest  eigenvalue  at  the  right  of  X  and  the  largest 

eigenvalue  at  the  left  of  X  respectively,  we  take  3.0  and  2.0.  The  initial  vector 

v^  is  generated  randonly.  After  each  outer  iteration  uaing  a  polynoaial  of  degree 

40  a  projection  step  using  the  last  six  vectors  is  taken.  The  Ritz  eigenvalues 

are  coaputed  and  coapared  with  those  obtained  at  the  previous  projection  step  as 

described  in  section  6.  When  either  of  the  eigenvalues  corresponding  to  X*.  X,  or 
-  -3 

X  has  converged  to  a  relative  accuracy  of  10  then  b,  p,  or  c  is  replaced 


accordingly,  when  the  three  eigenvalues  have  converged  the  approximate 
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eigenvector  is  replaced  by  the  Bitz  vector.  Figure  7-4  shows  480  steps  of  such  a 
process  (curve  1) .  Note  the  draaatie  drop  of  the  residual  norm  after  the 
projection  process.  Clearly  such  an  improvement  is  not  only  due  to  the  fact  that 
we  now  use  better  paraaeters  b  and  c,  but  also  to  the  fact  that  the  current 
approziaate  eigenvector  is  now  purified  from  the  coaponents  in  the  directions  of 
the  undesired  eigenvectors. 

As  a  coaparison  the  residual  noras  are  plotted  for  the  case  where  the  exact 
paraaeters  X  -  2.436872,  X  -  2.471196,  and  X+  “  2.604246  are  uaed.  In  this  case 
we  iterated  with  a  polynoaial  of  degree  100  (Iteration  with  a  polynoaial  of  degree 
40  gave  auch  slower  convergence)  .  Note  that  the  process  with  orthogonal 
projection  is  quite  successful  here  since  it  does  a  better  job  than  the  algoritha 
using  the  optiaal  paraaeters  and  a  reasonably  high  degree  polynoaial. 

8.  Conclusion. 

The  nuaericsl  experiments  suggest  that  the  use  of  orthogonal  polynomials  for 
solving  indefinite  linear  systeas  as  described  in  this  paper  can  be  effective 
especially  if  one  or  aore  of  the  following  conditions  are  net: 

-  The  systea  to  solve  is  not  too  badly  conditioned. 

-  A  aoderate  accuracy  is  required 

-  The  operations  y  ■  A  x  arc  vary  cheap 

-  Several  systeas  with  the  seas  aatrix  A  are  to  solved.  In  that  case  the 

paraaeters  are  estiaated  only  once  . 

The  Generalized  Chebyshev  Iteration  is  a  stable  process  and  relatively  high 
degree  polynoaials  can  be  used  without  difficulty  to  achieve  a  better  performance. 
In  order  to  estiaate  the  optimal  paraaeters  we  resort  to  a  projection  procedure 
which  incidentally  oan  also  be  used  to  improve  the  current  approziaate  solution  by 
reaoving  the  undesired  coaponents. 
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Several  problems  both  theoretical  and  practical  remain  to  be  solved  and  more 
numerical  experience  is  needed  before  this  class  of  techniques  mill  become 
reliable. 

It  is  also  hoped  that  besides  their  use  in  Numerical  Linear  Algebra,  the 
polynomials  introduced  in  this  paper  will  find  applications  in  other  fields  of 
numerical  analysis  such  as  in  approximation  of  functions. 
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